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Abstract 

Using a new parallel computing technique, we have run the largest cosmic string simulations 
ever performed. Our results confirm the existence of a long transient period where a non-scaling 
distribution of small loops is produced at lengths depending on the initial correlation scale. As time 
passes, this initial population gives way to the true scaling regime, where loops of size approximately 
equal to one- twentieth the horizon distance become a significant component. We observe similar 
behavior in matter and radiation eras, as well as in flat space. In the matter era, the scaling 
population of large loops becomes the dominant component; we expect this to eventually happen 
in the other eras as well. 



I. INTRODUCTION 



The formation and evolution of field theoretic cosmic strings has been extensively studied 
for the past thirty years (See [l| and references therein). The idea jij that superstrings 
could be stretched to cosmological scales by the expansion of the universe has recently [sHsl 
revived the interest in cosmic string networks. The existence of a cosmological network of 
strings may yield observational signatures that can be detected with current or planned 
experiments, giving us the opportunity to probe new physics at tremendous energy scales. 
This has motivated the study of many aspects of cosmic strings. 

Early work on cosmic strings focused on the simple and generic observational predictions 
arising from their kinematic gravitational effects, namely searching for their imprint on 
the cosmic microwave background via the Kaiser- Stebbins effect or by looking for the 
telltale identical pair of images of an astrophysical object being lensed by a cosmic string 
These effects are enhanced by increasing the energy scale of the string, and can therefore 



be used to place an upper bound on the dimensionless string tension G/x |8l4l0l|. Another 
observational signature predicted by the evolution of a cosmic string network is the existence 
of a stochastic background of gravitational waves [llj emitted by the oscillating string loops 
which continually break off the network. This is an important effect since it allows the loops 
to shrink and decay, preventing them from becoming a dominant contribution to the energy 
budget of our observable universe. Particular features of the string evolution, such as cusps 
and kinks, create a focused burst rather than stochastic gravitational wave emission, and 
so allow much lighter strings to produce a detectable intermittent signal [l2|. Additionally, 
cosmic strings can also produce other forms of radiation, either due to the partial annihilation 



of the string itself in regions of high curvature |l3l-ll6l| , or due to the their coupling to some 



other non gr av it at ional degree of freedom such as an axion [17[ , a dilaton [18| , or some other 



light field [19[. Radiation of this type has been studied in connection with cosmic ray physics. 



where we can use the observational bounds on these fluxes to place limits on the parameters 



of the cosmic string models [20 



However, it is clearly necessary to understand the statistical properties of the cosmic 
string network in detail before we can obtain robust predictions of observational signatures. 
Characterizing the essential properties of the network throughout its evolution has been 
approached in several different ways. Early work on this subject studied the properties of the 
network with analytical methods 21-2^|. The central thesis of this work was the approach 



of the network to a scaling solution whereby the strings contribute a constant fraction of 
the energy density of the universe. The existence of a scaling solution is paramount to the 
viability of cosmic string models, since the network would otherwise become the dominant 
fraction of energy in the universe, in clear contradiction with observations. 

After the initial impetus from the analytical description of the scaling solution, people 
turned to numerical simulation as a technique to compute the relevant parameters of the 
string network. Several groups independently developed codes to evolve a network of Nambu- 



Goto equations in an expanding universe |24|-|28| . Perhaps the most interesting conclusion 



of all these papers was that there does exist a scaling regime for long ("infinite") strings, 
where the average distance between the strings, d{t), and the coherence length, C,{t), both 
scale in proportion to the particle horizon distance, given by dh{t) = arj, where rj is the 
conformal time and a the scale factor. 

d{t) ~ ^{t) ~ 4(t) ~ t. (1) 
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In fact, although the codes were significantly different, they appeared to be close to a 
quantitative agreement on the parameters of the network. 



Further analytical work has been done by several different groups 29l-l3l|. but the ne- 
cessity of large numerical simulations in understanding the properties of the string network 
over a large range of scales has always remained. The constant improvement of computers 
and the introduction of new algorithms allowed for a second generation of numerical simula- 
tions 32-36] of much larger size. Ref. j33[ studied the loop number density and loop energy 
distribution and found an approach to scaling values of these quantities for loop sizes above 
a few thousands of the horizon size. Then Refs. HUel studied the loop production function 
and found an approach to scaling in that function, but their results were not compatible 
with the loop distribution found by (33| . 



Refs. |14I | simulated cosmic string loop networks using lattice field theory, but their 

results did not agree with the Nambu-Goto simulations of other groups. Instead, they found 
very significant emission directly from long strings. Attempts by other groups to reproduce 
this effect in field theory simulations 15|, |l6| did not succeed. 



One of the realizations which emerged from all these efforts was that there are two 
different time scales associated with the approach of a network to scaling. The first time 
scale parametrizes the way long strings in the network reach a steady state with the correct 
macroscopic properties of the scaling solution, and a second, much longer time scale is 
associated with the small scale structure of the network. In other words, the component of 
the network produced in loops, as well as the short wavelength spectrum of perturbations on 
long strings seem to approach scaling at a slower rate. This makes their study much more 
difficult in numerical simulations, where one is limited by computational resources. This 
has motivated us to develop a code based on new parallel simulation techniques which 
allows us to run simulations with dynamic range an order of magnitude larger than those 
previously reported in the literature. The results of these simulations will be described in a 
series of papers. Here we focus on the loop production function. 

With a larger dynamic range, it becomes possible to clearly distinguish behavior at two 
relevant length scales, namely the "smallest scales" (i.e., the simulation resolution, the 
gravitational backreaction scale, or in our case the scale of initial conditions), and the scale 
set by the particle horizon. There is an ongoing controversy as to which scale is important 
in determining the typical size at which loops are produced. Rou ghly speaking, the claims 



are that all loops are produced at scales set by the particle horizon 21, 22||, that all loops are 

or that some mixture of both are produced [3l| 



produced at the smallest scales |26|-|28 



38 



Our results confirm those of Olum and Vanchurin [36|, and our new techniques enable us 
to study loop production in much more detail and to much later times. We see a significant 
fraction of loop production in a broad, scaling peak reaching downward from about one 
twentieth of the horizon scale, along with an initially dominant but decreasing population 
at the smallest scales. This is true in both the matter and radiation eras, as well as in fiat 
space. In the matter era, the horizon-scale population becomes dominant at late times. We 
expect this to happen as well in the radiation and flat-space cases, but we are not currently 
able to run long enough simulations to see whether such regimes develop. There appear to 
be two different peaks in the scaling spectrum of loops at significant fractions of the horizon 
size. 
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II. SIMULATION TECHNIQUES 



The first step of any cosmic string simulation is to establisli tlie initial distribution of 



strings. We have used the procedure of Vachaspati and Vilenkin [39j to generate the initial 
configuration of the network. This procedure constructs a string from the square faces 
that it must pass through. For each such face, we choose a random position on that face 
equally distributed in a square filling the central 20% of the face in each direction, and we 
connect each position to the next by a straight string segment. For each such segment, we 
then choose a velocity v in a random direction perpendicular to the string with v'^ evenly 
distributed between and 0.25. These random perturbations mitigate the lattice nature of 
the initial conditions. 

The algorithm for the subsequent evolution of Nambu-Goto strings is based on that pre- 



sented in [32[, where the strings are described by a collection of straight segments linked 
together to form kinky long strings and loops. These piecewise linear strings can be sim- 
ulated exactly (up to computational arithmetic error) in flat space [32]. In an expanding 
background one can write the equations of motion for the string in conformal coordinates so 
that the causal structure of Minkowski space is preserved. The resulting equations cannot 
be analytically solved, but they can be well approximated by expanding the solution to first 
order in the product of if, the expansion rate of the universe, and the segment size being 
simulated [36|. This approximation improves as the simulation progresses since the Hubble 
distance grows relative to the size of the segments. 

The other important ingredient in this algorithm is the intercommutation of string when 
two pieces of the string worldsheet intersect. It is believed that most gauge theory cosmic 
strings would intercommute with a probability P of essentially unity. But this may not be 
true for fundamental strings in a 4-dimensional compactification, where the probability of 
such a process is suppressed. In the present paper we will only explore the case P = 1. 
Other cases will be studied in subsequent publications. 

The evolution of the string network continually produces loops of all sizes. It is useful 
to distinguish between the sub-horizon "loops" and the super-horizon "long strings" which 
make up the rest of the network. In our periodic space, all simulated strings are technically 
loops, and so we classify them by defining loops to be those strings which will never intersect 
any other string, including themselves. Due to the nature of our evolution algorithm it 
would be computationally intractable to evolve tiny loops, so we must remove them at some 
point. Loops whose physical size is much smaller than the interstring distance are likely to 
evolve without intersecting with any other string. In fiat space we expect small loops to 
fragment until eventually all the string energy in the loop ends up in non-self-intersecting 
trajectories. It is then safe to log and remove those loops from the simulation, since they 
would not have any effect on the network properties. In an expanding universe the situation 
is complicated by the fact that loops are affected by the cosmological friction, so a non- 
self-intersecting loop could, in principle, be destabilized and chop itself up, triggering a new 
stage of fragmentation. We expect this effect to be negligible for any loop of size significantly 
smaller than the Hubble distance. We have checked that this is not an important effect by 
letting non-self-intersecting loops oscillate a number of times before being removed. The 
results do not show any significant deviation in any of these cases. 

One of the main goals of this project is to extend the dynamic range of prior simulations 
performed using similar algorithms. Due to the fact that the simulations exist in a periodic 
space, the dynamic range is limited by the time it takes information to propagate across 
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the box. A possible solution to this problem was implemented in 32, 35[ by periodically 
doubling the size of the box of the simulation. Here we have chosen an alternative way to 
extend the running time, by splitting the box into different regions that are simulated by 
many different computers, in other words using a parallel code. This allows us to faithfully 
simulate cosmic strings in a much bigger box, and consequently to run to much later times. 
We developed, to this end, a parallel Nambu-Goto string simulation using the algorithm of 
32 . 35 , 36 1 and reported on preliminary results in j4o[. This parallel method proved to be 
indeed superior to the single processor codes, but it also had some inefficiencies associated 
with the very different work loads that different processors encounter in a simulation of this 
kind. At late times, some regions of the simulation volume will have little or no string, while 
others will have a great deal. In such a case, the power of parallel computing is greatly 
diminished, since all the other processors must wait for the one with the most simulation 
work to do. Faced with that difficulty we have developed a new parallelization technique, 
whereby the spacetime volume is divided into small 4- volumes whose boundaries are lightlike, 
i.e. they have only initial and final surfaces. This permits the individual 4-volumes to run 
independently of one another, imposing only that each individual job must wait for those 
who provide data for its initial surfaces. This new method is particularly useful for our type 
of simulation, and it has allowed us to increase the dynamic range by almost an order of 
magnitude with respect to all previous cosmic string simulations. A detailed account of this 



new technique and its comparison with conventional parallelization methods is given in [37 

Our technique naturally gives rise to a periodic simulation volume in the form of a rhombic 
dodecahedron with opposite faces identified. A point and its images form a face-centered 
cubic lattice. We will refer to the constant comoving distance between a point and its closest 
images as the size of the simulation volume, L. The comoving volume being simulated is 
thus V = L^/\/2. The face-centered cubic lattice is close-packed, so the ratio of V to Li^ 
is the smallest possible. As shown in j36|, the rate of information propagation through the 
simulation volume is about half the speed of light, so in a box of size L, we can run a 
simulation whose duration in conformal time is also L. 

We have run simulations in fiat space of size 2000 for a conformal time duration0 2000 
in units of the initial cell size of the Vachaspati-Vilenkin algorithm. Such a simulation 
contains at maximum about 14 billion linear pieces of string simultaneously, and in the 
course of evolution creates about 1 trillion segments and produces about 10 billion loops of 
string. If one counts the largest total amount of data existing at one moment of simulation 
time, we believe this ranks among the largest simulations of any kind ever performed. 

In the expanding universe we are not able to simulate as large a volume. The reason 
is that as the universe expands, the comoving length of linear pieces of string decreases 
(i.e., they are not purely stretched with the expansion). Thus, while the density of string 
in expanding cases is not very different from the fiat case at the same time, the number of 
pieces of string to be simulated is much larger, and so the simulation effort much greater. 

In the case of the radiation era, this has limited us so far to simulations of size 1500. In 
the matter era, the expansion is much more rapid and the problem more severe; the largest 
simulations we have so far run had size 500. 

Once one has decided to create strings with a Vachaspati-Vilenkin cell size of 1, one must 
choose a conformal time rji for the start of the simulation. This time determines the horizon 
distance, which will be used to define scaling, and in the expanding universe it controls the 
Hubble constant, which determines the rate of redshifting and stretching. We choose r]i so 

^ We work throughout in units where the speed of hght is 1. 
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that the initial conditions correspond as closely as possible to the conditions that one would 
expect (based on later stages in the simulation) in a scaling solution at that time. 

In the matter era, we found the best initial conformal time to be 4.5. We can thus 
run until conformal time 504.5, without contamination from the periodicity of the volume. 
Our dynamic range, defined as the ratio of final to initial conformal time, is thus 112. In 
the radiation era, we chose initial time 6. We can thus run until conformal time 1506, for 
dynamic range 251. In fiat space, we chose initial time 4.0 and ran until time 2004, for a 
dynamic range of 501. 

These dynamic ranges are much larger than those of any previous simulation. The largest 
simulation before this was done by our group [SG*] and claimed dynamic range 60 in the matter 
era and 120 in the radiation era. However, that claim was based on an artificial setting of 
the initial clock to 2.0 and 1.0 respectively. The initial conditions in those simulations were 
more appropriate for initial times of order 4.5 and 6.0, giving effective dynamic ranges 27 
and 20. Ref. [33[ used dynamic range 8 in matter and 17 in radiation, and all other previous 
simulations were smaller. A snapshot of a flat-space simulation is shown in Fig. [TJ 



III. SCALING NETWORK DYNAMICS 

Suppose that the cosmic string network does evolve into a scaling regime. Then we hope 
that simulations will show an approach to this regime. What should we look for? 



A. String density 

First, let us consider scaling of the long string energy density, which is the easiest to find. 
Let d denote the average interstring distance of long strings, defined by 



d=^—, (2) 




where /i is the tension of the string and poo is the energy density of the long string network. In 
a scaling regime we expect this distance to be proportional to the horizon distance dh = arj, 
so that 

is a constant. 



B. Loops 

Now let us turn to the distribution of loops. We characterize loops by their energy fil, and 
we call / the length. Some groups have concentrated on the the loop density distribution, 
which we will denote here aqj 

/, , dN 



^ But note that 



35 



36| used n(l,t) to denote the production density, which we caU f{l,t) here. 
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FIG. 1. A picture of the string network at time 500 in a flat-space simulation of size 500. Long 
strings (here, any loop longer than the horizon size) are shown in light gray and loops are colored 
(shaded) according to their length, defined as total energy divided by /i. Non-self-intersecting 
loops smaller than length 10 have been removed. The edges of the simulation volume are shown as 
black straight lines. The rhombic dodecahedron is seen from one of the order-4 vertices, so that its 
projection looks square. The side of the square is the horizon distance, 500 units. The largest loop, 
shown in red (dark gray) on the far right above the centerline has length about 148. It appears 
much smaller because it is wrapped into a closed loop, because it has depth that cannot be seen, 
because it is wiggly, and because its length includes its kinetic energy. Loops of this size are rarely 
seen in the loop production function, so this loop will probably fragment or rejoin the long string 
network. 



so that n{l, t)dl is number of loops in the network at time t whose lengths lie between / and 
/ + dl. Others concentrate on the density of loop production, which we will denote 

so that f{l,t)dtdl is the number of loops produced per unit volume with I E [1,1 + dl] at 
times t G [t,t + dt]. 

To put these functions in terms of scaling variables, we trade the loop length / for a length 
in scaling units, x = l/dh- Similarly, we will use a scaling length interval dx and a scaling 
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spatial volume df^ or spacetime volume li^. Thus we define the dimensionless functions 

n{x) = 4n{l,t) (6) 

and 

f{x) = dlf{l,t). (7) 

We write the left-hand sides as functions of x, since in a scaling regime these functions will 
depend only on x, not on time. 



C. Energy balance 

The production of loops and scaling of long strings are related. If we consider a long 
piece of string whose length is I, the expansion of the universe will stretch and redshift it 
according to the relation [22| 

|=i7/(l-2(.L» (8) 

where (f^) is the average squared velocity of the long strings. The energy density in long 
strings then obeys the Boltzmann equation 



dpc 



dt 



POO 

-2H{l + {vl))p^-fi lf{l,t)dl. (9) 

Jo 



Now consider a cosmology with scale factor a ~ t'', so that u is equal to in fiat spacetime, 
1/2 for radiation domination, and 2/3 for matter domination. The horizon distance is then 
t/{l — u). In a scaling regime poo = p/{l'^d\), so we can rewrite Eq. ^ as 



mM)c^/ = ^^i-^(0)- (10) 



In scaling variables, Eq. ( jlOl) becomes 



xf{x) rfx = (l - Y^ivl)^ ■ (11) 

A consequence of Eq. (ITT]) is that any scaling loop production function f{x) must be 
normalizable, in the sense that the left-hand side of Eq. (ITT]) converges. 

Below, we will measure 7 and (f^) in our simulations and use them to predict the left- 
hand side of Eq. fill I) , as a check on the consistency of our results. 



D. Integrated loop production 

The loop density distribution arises from all prior production of loops. We will ignore 
gravitational damping, which is not included in our simulation. We also, at first, ignore the 
decrease in the energy of loops from redshifting of their center-of-mass velocities. In this 
approximation, sub horizon loops retain their physical length /, and the only change in the 
distribution of loops in a comoving volume is due to the production of new loops from the 
long string network. Thus 
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j^[a%t)n{l,t)]=a'{t)f{l,t). (12) 
Accordingly we can find n{l,t) from f{l,t) by integration, 

n{l, t) = J(f)l a\t')f{l, t') dt' = £ (^0 /(/, t') dt', (13) 

since a{t) ~ t". 

We would like to write this in terms of the scaling functions n{x) = d'^{t)n{l,t) and 
f{x) = d\{t)f{l,t). Let us change the integration variable from t', the creation time of the 
loop, to x' = l/dh{t'), the scaling length that it had when it formed. We similarly exchange 
t' /t for x/x' to get 

/■oo 

n{x) = (1 - u)x^''-^ / x"^-^'f{x') dx'. (14) 



Since v is at most 2/3, 3 — 3i/ is at least 1. Since f{x) must be normalizable (in the sense 
of Eq. (ITT]) ), the integral on the right hand side of Eq. (IT^ converges even if x is taken to 
0. Thus for loops sufficiently below the horizon size {x ^ 1), Eq. ( IT^ is insensitive to the 
lower limit, and we get a power-law prediction [T] for the loop density distribution, 

n{x) = u^x^''-^ (15) 

with 

poo 

uj^ = {l-iy) / x^-^^'fix) dx. (16) 



JO 

Now let us consider the effect of the center-of-mass velocity of emitted loops on Eq. ([1] 
We will continue to neglect gravitational damping. As we will show below, typical loop 
velocities increase with decreasing x, and tiny loops are usually emitted with quite large 
boosts. For simplicity, we will assume that the dominant contribution to Eq. f|T^ comes 
from a single Xi corresponding to a typical loop center-of-mass velocity Vi. For x near xi, 
the effects are somewhat complicated [4l|, and we will not attempt to model them here, but 
for X <^ Xi matters are simpler. In this case, the most important loops have been redshifted 
essentially to rest, and thus their energy has decreased by a factor 71 = (1 — u^) 

Thus we can update our calculation by taking loops with length / to have been formed 
with length I' = 71L Equation (|T3|) becomes 

—f{l',t')dt'. (17) 

The factor of 71 arises because n{l, t) is the loop density per unit /, while /(/', t') is the loop 
production density per unit /' and dV jdl = 71. 

The relationship between x' and t' is now x' = l'/dh{t') = •yil/dhit') so x/x' = t'/(7it), 
and we find once again that for x ^ xi we have Eq. f|T5l) . but with 

PCX) 

uj, = {l- z/)7^"^ / x'^-'^^fix) dx. (18) 
Jo 

Equation f|T5|) tells us that the scaling loop density distribution n{x) has a universal form 
(for small x) that does not depend on the shape of f{x), except in an overall factor. The 
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only way that a simulation could find a function n{l, t) whose / ^ t behavior does not have 
the form (i^^n(x) with n{x) given by Eq. (fTS!) is for the loop production in that simulation 
to not be in a scaling regime. 

This conclusion applies even though real cosmic strings, unlike those in usual simula- 
tions, would be affected by gravitational damping. It does not make sense to claim that 
simulations indicate a nonnormalizable f{x), but that f{x) will be cut off for low x by gravi- 
tational damping. If simulations indicate a scaling function f{x), then if one did a very long 
simulation, one should observe f{x) approaching that scaling form. In order to approach a 
nonnormalizable f{x), the fraction of the energy of the string network emitted into loops in 
a Hubble time in the simulation would have to grow without bound. Such a situation could 
never be found in a simulation that conserved energy, so such a result does not make sense. 

The same argument applies to analytic claims about loop production functions. It is 
not reasonable to predict a nonnormalizable f{x) on the basis of analytic reasoning that 
ignores gravitational back reaction, and then claim that such a thing is made acceptable by 
the inclusion of gravity. What could such analysis predict about strings with infinitesimal 
Gfi, so that gravity could be ignored? Some consistent answer should emerge, which means 
that the nonnormalizable growth of f{x) would have to be cut off at some small x. One 
must then analyze whether this cutoff is at a smaller or larger scale than that produced by 
gravitational damping, which would depend on the value of Gfi under consideration. 



IV. RESULTS: LONG STRING DENSITY 

We now present results from our simulations, beginning with the scaling of the long string 
density. Long strings are commonly defined as those whose energy is above some threshold. 
Here, however, we distinguish the long strings from the loops based on interactions. If a 
loop undergoes 3 oscillation^ without intersecting itself or another string, we retroactively 
consider it to be a loop beginning at the time when it was formed. Everything else belongs 
to the long string network. 

With this definition, we compute the parameter 7 = d/dh, with d given by Eq. ([2]). For 
a scaling network we expect 7 to achieve a constant value. In Fig. [2] we plot 7 as a function 
of conformal time 77. The top group of lines shows 13 runs of size 500 in the matter era, the 
middle group 6 runs of size 1500 in the radiation era, and the bottom group 4 runs of size 
2000 in flat spacetime. In all cases, the scaling of long strings is well established early in the 
simulation. We note that in the fiat case, the scaling results entirely from loop production, 
as there is no Hubble friction present. 

Nevertheless, as seen in this plot against logarithmic time, the interstring distance is 
slowly changing, even at end of the simulations, so the final value cannot be precisely 
determined. At late times, there is also quite a lot of noise because of fluctuations in the 
actual long loops produced. 

We average the interstring distance over all runs in each era at the end of the simulation, 
and show the results in Table HI in comparison with previous simulations. In general we see 



good agreement with other recent results. The exception is the field theory simulations [10 



which do not seem to agree with any Nambu-Goto simulation, and the flat space simulations 



of Ref . [3^ , which gave an interstring distance twice as large as ours. 



^ This is the minimum number of oscillations before our algorithm can confirm that the trajectory is non- 
self-intersecting. Changing this to a larger number does not significantly affect the results. 



9 



0.3 




10 1 00 1 000 

FIG. 2. The ratio of the interstring distance to the horizon size. The top group of hnes is for the 
matter era, the middle group the radiation era, and the bottom group flat spacetime. 



First author &: Ref. 


Flat 


Radiation 


Matter 


Albrecht [25] 




0.07 


0.12 


Bennett [27] 




0.14 


0.18 


Allen [28] 




0.13 




Vanchurin [32] 


0.096 






Ringeval [33] 




0.162 


0.188 


Martins [34] 


0.23 


0.13 


0.20 


Bevis [10] 




0.255 


0.285 


This paper 


0.12 


0.15 


0.17 



TABLE I. Values of ^, the ratio of the interstring distance to the horizon size, from present and past 
simulations. Ref. [lOj] used lattice field theory and included all string in y^Qo . All other simulations 
used Nambu-Goto, and all but the present paper included loops larger than some scaling threshold 
in Poo- Here we include all loops which will fragment (or rejoin) in p^o- This is the correct 
definition for the calculation of energy conservation in WDI below. Including only loops larger 
than the horizon size would increase 7 by no more than 3%. 

The agreement with other simulations whose dynamic range was significantly smaller 
than ours supports the conclusion that the long string component of the network finds its 
way to a scaling solution on a relatively short time scale, even if the other properties of the 
network have not relaxed to their scaling solutions yet. We will comment on this effect again 
when we discuss the scaling of loops. 
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V. RESULTS: LOOP PRODUCTION 



In this section we study the loop production function /(x) and the loop distribution 
function n{x), as defined above. We consider a loop to be produced when it first enters a 
non-intersecting trajectory. Loops that fragment or that rejoin the long string network are 
not counted. 

The distribution f{x) dx gives the spectrum of loop production, which is not necessarily 
integrable. Instead we seek the distribution of power going into loops xf{x) dx. Due to the 
wide range of loop sizes covered by our simulations, it is illuminating to plot this power 
density in logarithmic units: x'^ f{x) dlnx. In order to interpret the area under the curves 
as the total power produced in loops, we plot x"^ f{x) on linear-log axes in the figures below. 

A. Matter era 

We have performed thirteen simulations on a comoving box of 500 initial correlation 
lengths each, in a matter era background. The initial conformal time was set to rji = 4.5, 
and we ran until rif = 504.5. Following the algorithm described above, we removed the 
non self-intersecting loops from the simulation, keeping record of their size and time of 
production. Using this information, we can reconstruct the loop distribution n{x) and the 
loop production function f[x) as time progresses to determine if either approaches a scaling 
regime. 

We first turn our attention to the loop production function. We plot in Fig. [3] the data 
for f{x). The error bars show one standard deviation from the thirteen successive runs. The 
form of the curves in Fig. [3] is similar to that of [sHi , but because our simulation is larger 
than that of [36] we can follow the evolution of f{x) to much later times. 

Notice these curves have a slightly increasing area toward later times. This suggests the 
interstring distance d began at a value larger than its scaling value, and the initially lower 
power going into loops allowed it to decrease toward the scaling value. This is confirmed 
in Fig. [21 Why didn't we choose our initial time so that the initial •y = d/dh was equal to 
the scaling value? The reason is that we found it more important to have the the height of 
the small loop production peak start out unchanging in time, so that we could be sure that 
the decrease at later times (shown in Fig. [3]) was not an artifact of the initial conditions. 
One cannot make both choices simultaneously because the structures on strings in the initial 
conditions are not very close to the scaling regime. 

We can get a better look at the relative shapes of the curves in Fig. [3] by normalizing 
each curve to unity. This is done in Fig. HJ which now gives the relative flow of power into 
loops of different sizes. We see that the normalized f{x) appears to converge more rapidly 
to a final scaling form. 

We can identify three different peaks in the loop production. The first is sharply peaked 
at x ~ 0.05. The second is centered approximately an order of magnitude smaller in x, and 
is much wider. The third peak is clearly moving toward smaller x, representing the transient 
(non-scaling) population of loops seeded by the initial conditions. The non-scaling peak is 
decreasing in area relative to the scaling portion of the distribution, and is subdominant 
after a conformal time of order 100. At the farthest reach of our simulation, the non-scaling 
power represents of order 25% of the total loop power, and we expect it to continue to 
decline in larger simulations. 

It is clear from the time dependence of f{x) that the small scale structure and loop 
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FIG. 3. The spectrum of loop production power f{x) in the matter era. 



production has not completely reached scaling, but we consider the late-time behavior of 
our data to be evidence that a time-independent loop production function exists without 
the aid of an additional smoothing mechanism such as gravitational radiation. While it 
is possible that the non-scaling third peak will continue to evolve toward smaller x while 
maintaining a constant (but certainly subdominant) area, we do not expect this to be the 
case. 

The other two peaks appear to scale: they do not move toward smaller x at later times. 
We do not know why there are two such peaks, but one speculation is that one peak rep- 
resents loops which resulted by chance from two or more spacelike-separated intercommu- 
tations, while the other represents those that were formed from features on a long string 
by a single reconnection. Another possibility is that the peak with smaller x represents the 
results of fragmentation after loops are formed, while the peak with larger x represents those 
which happen to form in non-self-intersecting trajectories. 

To study the center-of-mass velocities of loops, we define p = v/\/l — v"^, the loop mo- 
mentum per unit rest mass, and let /(x, p) dx dp be the scaling production rate of loops with 
X G [x, X -|- dx] and p G [p,p + dp] . To show x/(x, p) dx dp we plot x'^pf{x, p) on logarithmic 
axes in Fig. [51 As expected, small loops at late times are mostly formed with ultrarelativistic 
speeds. 

We now plot in Fig. [6] the number density distribution n(x) dx, which is represented by 
xn(x) on logarithmic axes. We see that the solution is in good agreement with the results 
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FIG. 4. Normalized f{x) during the matter era. While perhaps 70% of the power goes into 
initial-conditions loops at conformal time 50.0, less than 25% does so by conformal time 500.0. 

of Eq. f lTSj) . which for the matter era predict a power law behavior of the form x^^. This 
is different from the result obtained in j33|, where the power law was found to be of the 
form x"^ "^. As we showed in §111 D I above, the scaling form of xn{x) could not go as 
throughout the range of x where loops are produced. Either there must be significant loop 
production in the simulations of [ssj at scales below where the x~^''^ fit applies or their loop 
distribution must not have reached its final scaling form. 



B. Radiation era 

Fig. [7] shows less dramatic, but similar, scaling behavior in the radiation era, where we 
performed six simulations of size L = 1500 comoving initial correlation lengths. The results 



are again qualitatively similar to those of (36| . Notice that the total power going into loops 



is much higher for radiation domination than for matter. 

The approach to loop scaling is much slower (in conformal time) than occurs in our 
simulations of the matter era. Nevertheless, it appears that a large but shrinking population 
of transient loops at the initial conditions scale is making way for a scaling population sharply 
peaked at a; = 0.05, with a nearly flat distribution extending over many orders of magnitude, 
perhaps peaked an order of magnitude smaller in x. 

We can understand the slower approach to scaling as follows. During radiation domi- 
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FIG. 5. Contour plot showing the distribution of loop power in the matter era, for conformal times 
from 283.1 to 504.5. Here p = vj^JX — . 

nation, the increase of scale factor with conformal time is much lower, and so strings are 
being stretched and thus smoothed more slowly, causing small-scale structure to persist for 
longer. Furthermore, the redshifting of string energy is less efficient than during the matter 
era, and therefore the network will need to make more use of loop production in order to 
reach a scaling solution, so the overall production rate is much higher. 

We cannot predict the scaling form of j[x) as well as we could for the matter era, since 
a dominant fraction of power is in the transient regime. But the late-time radiation era 
power looks similar to that in the matter era at a much earlier time. We expect that the 
radiation era evolution should be similar to that in the matter era, with the scaling features 
eventually dominating, and /(x) eventually approaching a time- independent form, but the 
approach to scaling will last much longer. 

The velocity dependence of the loop production function can be seen in Fig. |H1 
We show in Fig. |9] our dataset n{x) binned logarithmically from the simulation. The 
data obtained this way agrees quite well with the analytical results presented earlier for the 
radiation era, where xn{x) = ujix~^^'^. This should not be misinterpreted as true scaling, 
since we know that a dominant portion of the power in loops is in a transient, not a scaling 
population. 

Roughly speaking, there are three possibilities which can occur at very late times (as- 
suming 7r is static). The non-scaling peak could smear itself out as it continues to move 
toward smaller x while creating an extremely broad plateau. It could conceivably move off 
the plateau and continue indefinitely toward smaller x while maintaining a constant area. 
It could also decline while the peak at a; = 0.05 grows, and not leave a tremendous plateau. 
We conjecture the last to be the case. It should be pointed out that the first two cases 



14 




FIG. 6. The number density distribution of loops xn{x) during the matter era. 



are essentially identical in terms of the prediction for ui , but that third case will cause a 
relative increase by up to an order of magnitude. It is this uncertainty which prevents us 
from claiming that the loop density shown in Fig. M is in the final scaling regime. This 
illustrates why the spectrum of loop power f{x) is a more precise indicator of scaling than 
n{x), the spectrum of existing loops. 



C. Flat space 

We present in Fig. [10] the results of our four simulations performed in a box of size 2000. 
Because redshifting is nonexistent, all smoothing of the strings comes from loop production. 
This massive loop production clearly smooths the strings enough for a scaling peak at x ~ 0.1 
to appear, although it remains small even until conformal time 2000. 

The results in Fig. [10] are qualitatively similar to those of [35], but the ratio of the 
scaling to the non-scaling peak height in [35i] was si gnifi cantly larger. There are several 



differences between the present simulation and that of j35j, which make a direct comparison 



difficult. Ref. [35| achieved box size 800 by a technique of successive doublings of the box size 
3^ . whereas we simulate a box of size 2000 directly. The box-doubling technique required 
intercommutation probability P = 0.5, whereas we use P = 1. We run only for an interval 
of conformal time equal to the box size L, whereas [s^ ran significantly beyond conformal 
time 800. Finally, [351] used somewhat smoother initial conditions in which a point was 
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FIG. 7. The loop production power f{x) in the radiation era. 

interpolated between each pair of successive Vachaspati-Vilenkin faces, whereas we draw a 
straight hne between those points. Smoother initial conditions may have led to a lower peak 
at the initial-condition scale. Note also that jH, [SGj show x^f{x) on a logarithmic scale, 
whereas the vertical axis in Fig. [10] is linear. 

Perhaps the most important features of Fig. [10] are the turnaround of the transient peak 
and the existence of a scaling peak. After a time t] k, 150, the amplitude of the transient 
peak begins to decrease. Following the arguments presented above, we again expect that this 
decrease will continue, leading to a loop production function consisting of a time-independent 
distribution peaked at x ~ 0.1. 

D. Energy balance check 

Having measured values for 7 and (f^) directly from the simulations we can now use 
Eq. fill I) to predict the total power emitted by the long string network in the form of loops, 
^Prediction, and comparc it with the direct computation of this power from the integral of the 
loop production function, namely, 

poo 

^Simulation = / X f {x) dx. (19) 

Jo 

We show in Table [TT] the comparison of these quantities for the matter, radiation and flat 
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FIG. 8. Contour plot showing the distribution of loop power in the radiation era, for conformal 
times from 475 to 1006. Here p = vj^JX — v^. 
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^Prediction 


^Simulation 


Matter 


0.17 


0.35 


21 


19 


Radiation 


0.15 


0.40 


53 


51 


Flat 


0.12 


0.45 


139 


136 



TABLE II. Predictions from energy balance. 



eras. We see that the agreement between the predicted result and the numerically found 
value is within the statistical noise in all cases. We consider this is an important check for 
our simulations since it serves as a nontrivial test of our code and algorithms. 

VI. CONCLUSIONS 

We have developed a new parallel computing technique js^] that has allowed us to run 
the largest cosmic string simulations ever performed, reaching an increase in dynamic range 
of roughly an order of magnitude with respect to previous studies. The results of our 
simulations indicate a much slower approach to scaling for loops than for the long string 
contribution to the network. This clearly shows the need to push the dynamic range of 
the simulation to its largest possible extent in order to prevent contamination by transient 
states. 

How is it possible for the interstring distance to be quite close to the final scaling value 
at early times, while the loop production function is still very far from scaling? Surely 
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FIG. 9. The number density distribution of loops xn{x) during the radiation era. 



an important part (in flat space, the only part) of the loss of energy from the long string 
network is due to the production of loops, so if the loop production mechanism is not in a 
scaling regime, why does the energy loss rate scale? 



The explanation is presumably that described by Bouchet |42|. The loops are produced by 



a two-step process. The first step is the intercommutation of infinite strings, which maintains 
the interstring distance by the usual feedback mechanism in which an increase in the string 
network density leads to more intercommutations. The second step is the formation of loops 
triggered by or made up of the large kinks formed in the intercommutations. It controls 
how the removed string ends up in loops of different sizes. Thus the first step can be scaling 
long before the second. 



Our simulations confirm the results of [35|, |36| : There is an initially dominant transient 



feature in the loop production function which very slowly subsides to reveal a significant 
fraction of loops produced in a scaling regime. In the matter era, we find this population 
dominates over non-scaling loop production. Because this does not occur until extremely 
late times, we do not believe that prior claims of having found the final form of loop scaling 
in numerical simulations are correct. 

In light of our results, we expect that the final loop production function will eventually 
be dominated by a broad peak of loop sizes proportional to the horizon size but ranging 
downward from about a twentieth of that size. We achieve this to good accuracy only in the 
matter era, where approximately 75% of loop power has achieved scaling by a conformal time 
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FIG. 10. The loop production power f{x) in flat spacetime. 

of 500. This will clearly have an impact on the observational signatures of string networks. 
We will report on this in future pubhcations. 
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Appendix A: Infinitely difficult cosmic string loops 

We have identified some new features of Nambu-Goto cosmic strings which are impractical 
to faithfully simulate. One of these features can be called a "skipping stone." It occurs in 
particular cases when a network contains a loop which to good approximation is piecewise 
linear with five or six kinks, and which collides with a long straight segment of string. The 
collision may occur in such a way as to break off a loop again immediately afterward, one 
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with less energy but precisely the same shape. 

The problem is that the smaller loop will then repeat this process, skipping off the long 
string again and again, losing energy each time, but never changing shape. Just as a skipping 
stone leaves behind a geometric series of ripples, this loop will perform (in the Nambu-Goto 
approximation) an infinite number of intercommutations, leaving behind a geometric series 
of kinks on the string. 

Such a physical process represents a nightmare for a numerical simulation which has no 
minimum resolution, since each of these ripples will be recorded and evolved. 

To avoid the tremendous computational resources required to simulate these rare skipping 
stones all the way down to the minimum size set by the floating point resolution, we have 
intentionally failed to perform a certain, very small number of intercommutations. These 
avoided intercommutations are chosen to only include the collision of a disjoint loop with 
another string, and so never prevent the formation of a loop. Furthermore, the avoided 
intercommutation must occur within a very short distance of a very large number of kinks, 
i.e., only after the stone has skipped at least 50 times do we allow it to pass through the 
surface of the water. The fraction of intercommutations we ignore is less than 0.2%. 

Although these features occur with a varying degree of severity, and rarely plague small 
simulations, in a simulation of size 2000 they are virtually guaranteed to occur at least once 
with enough alignment so as to bog down one computer for literally days as it attempts to 
evolve through an amount of simulation 4-volumc which normally takes a few seconds. 

If cosmic strings exist, we expect that "infinite repetition" phenomena such as we discuss 
here occur in the real universe. Real cosmic strings would not have exactly straight segments, 
but this is of little consequence. If an approximation of the "skipping stone" phenomenon 
began, each successive iteration would involve shorter pieces of the original strings, which 
would thus be closer to straight. Of more importance is that the kinks separating the 
straight segments would not be infinitely sharp. Some kinks appear anew each cycle, but 
some remain from the original shape of the loop, and these would have been smoothed by 
gravitational back-reaction. Thus the gravitational damping scale would set a lower bound 
on the size of self-similar loops that could be produced. Since the number of breakings off 
and rejoinings grows only logarithmically with the ratio of the loop size to the curvature 
radius at the kink, and since the entire phenomenon is quite rare, we do not expect it to 
have any observable consequences. 
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